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1. Introduction 


The phenomenon of fatigue-crack-closure, first discovered experimentally 
by Elber [l,2], continues to be a subject of several recent experimental 
and analytical studies; see, for instance, Refs. [3,4]. Elber has also 
originally postulated that the crack-closure phenomenon is caused by residual 
plastic deformations remaining in the wake of the advancing crack- tip. 

Analytical models that lend theoretical support to the existence of 
the crack-closure phenomenon in fatigue crack growth, and provide some 
rationality for the adoption of an effective stress-intensity range, based 
on closure effects, for the correlation of fatigue crack growth rate, have 
also been proposed by Budiansky and Hutchinson [5]. 

As for a more general analysis of extending cracks under general block 
cyclic loading, to obtain crack-closure stresses, crack-opening stresses, 
details of crack-surface displacement, and residual stresses in the crack- 
tip region, etc., elastic-plastic finite element analyses were first per- 
formed by Newman and his colleagues [6,7,8]. Apart from these analyses, the 
authors are aware of similar attempts only by Ohji and his co-workers [9,10]. 
The studies in [6-10] considered the Mode I case only. Also, since the crack 
growth was simulated in [6-10] by shifting a finite element node (the current 
crack-tip) to an immediately adjacent node, and since constant strain tri- 
angle type finite elements were used to model the cracked structure, a 

very fine finite element mesh (with the smallest element often being of the 
-3 

order of 10 times the crack length) is necessary in the modeling procedures 
of [6-10]. This can be very expensive especially when cyclic loads of 
arbitrary spectrum are considered. 

One of the objectives of the present report is to present an alternate 
cost-efficient and accurate elastic-plastic finite element procedure to 
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analyse fatigue crack closure and its effects under general spectrum load- 
ing. Both Modes I and II type cycling loadings are considered. Also 
presented in this report are the results of an investigation, using the 
newly developed procedure, of various factors that cause crack growth 
acceleration or retardation and delay effects under high-to-low, low-to- 
high, single overload, and constant amplitude type cyclic loading in a 
Mode I situation. Further, the results of an investigation of a center- 
cracked panel under external pure shear (Mode II) cyclic loading, of con- 
stant amplitude, are reported. 

2. A Synopsis of the Present Analysis Procedure 

Before embarking on a presentation of the mathematical details, the 
salient features of the presently reported analysis procedure are given 
be low : 

(i) The present elastic-plastic finite element procedure accounts for 
arbitrary strain-hardening material behaviour; the mathematical description 
of incremental elastic-plastic flow consists of Huber-Mises-Hencky yield 
criterion, and a Prager-Ziegler type kinematic hardening rule which best 
models the Bauschinger effects. 

(ii) The well-known Hutchinson-Rice-Rosengren [11,12] type strain and 

stress singularities, for strain-hardening materials, are embedded in specially 

developed elements near the crack- tip. This eliminates the need for a very 

fine mesh near the tip. For instance, the crack- tip elements in the present 

procedure are of the order 10 ^ of the crack length, as compared to constant 

-3 -4 

strain triangles of the order of 10 to 10 times the crack length generally 
used in the procedures of [6-10]. A hybrid displacement finite element 
method [13,14] is used in developing these special elements. 
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(iii) The above special hybrid elements are of circular-sector shape, 
centered at the crack-tip; thus enabling one to model crack growth in any 
arbitrary direction, from the crack-axis, under general mixed mode cyclic 
loading. 

(iv) In the present procedure, crack-growth is simulated by: (a) 

translation of the core of circular-sector elements, with embedded H-R-R 
singularities, by an arbitrary amount in the desired direction, (b) reinter- 
polation of requisite data in the new finite element mesh , and (c) proportional 
relaxation of tractions in order to create a new crack surface. 

(v) In order to most accurately determine the crack-opening stress a 

and crack-closure stress cr , the displacements at nodes on the crack-axis 

cz 

before closure (and after opening), as well as the restraining force at the 
corresponding nodes after crack-closure (and before opening) are extrapolated 
against the load level. In all the cases studied, these two sets of extra- 
polated values for CT and cr^ were found to correlate excellently. 

(vi) A static-condensation procedure is employed wherein the plastic 
portion of the structure is isolated from the elastic; the stiffness of only 
the former part keeps changing whereas that of the latter remains fixed. 

This results in a considerable saving of the computational time. This and 
the invocation of the appropriate anti-symmetric and symmetric conditions 
for modes II and I problems, along with the use of special elements of a 
large size as described in (ii) above, reduce the analysis procedure to be 
feasible as a routine tool in design, if necessary. 

(vii) A study is made to arrive at a criterion for the stress level 

cr , at which fatigue crack growth occurs. In prior literature, this crack 
extension stress level was chosen arbitrarily. For instance, in [6-8] the 
crack is extended at the maximum applied stress in each cycle in a spectrum 
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loading, whereas in [9,10] the crack was extended at the applied stress 

level at which the restraining nodal force at the new crack- tip becomes 

zero. In the present study, for instance in a constant-smplitude (zero 

to tension) cyclic loading, it was found that O' and O ^ were very sensitive 

to the chosen o . In the present work, a criterion, <J = p(cr - cr ) 
ex ex max op 

+ cr where p is a constant of proportionality, is postulated; and p is 

obtained by calibration such that the calculated cr correlated with that 

op 

observed in experimental studies such as in [1,2,33. However, it is not 
claimed that this is a general criterion; it is recognized that it is dependent 
on material properties, and to an extent on the analysis procedure itself. 

The mathematical details of the above features are discussed in the following. 

3. Elastic-Plastic, Embedded Singularity, Incremental Finite Element 
Method, Based on Hybrid Displacement Model, For Analysing Cracked Structures 

To start with, we use circular-sector shaped "singularity" elements 
near the crack-tip, as shown in Fig. 1. In these singularity elements, 
a displacement field which corresponds to strain and stress singularities, 
for strain-hardening elastoplastic materials, of the well-known Hutchinson- 
Rice-Rosengren [ll,12] type, is assumed. The above singularity elements 
are surrounded by "regular" eight-node isoparametric quadrilateral elements, 
as shown in Fig. 2. Compatibility of displacements and continuity of trac- 
tions between these "regular" and "singular" elements is enforced through 
a Lagrangean Multiplier technique as shown below. 

The incremental analysis of the present elas to-plastic problem is 
based on a continuously updated Lagrangean (coordinate-system) formulation. 

In the following, for simplicity, we consider the formulation in the context 
of a general three-dimensional problem and present only the essential 
mathematical details. We consider a fixed cartesian coordinate system, and 
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consider to be the 'state' (viz. deformation, strain, and stress) of the 

structure before the addition of the Nth load- increment . Let x. be the 

1 

material coordinates of a point in the initial, undeformed and unstressed 

N 

configuration, C^. Let u^ be the displacement, measured in the fixed 

cartesian system, of a point in to the deformed state C^. The new 

coordinates of the material point in C„ then become X*? = x. + u^? . Let 

N xii 

N 

the (symmetric) Eulerean (true) stress tensor in be T„, measured per 

unit area in and in the metric of the fixed cartesian system. Let the 

body forces/unit volume in be and prescribed surface traction/unit 

area in C be T*!. 

N 1 

Let additional body forces Af\ (per unit volume) and additional 

surface tractions AT\ (per unit area) be applied during the movement of 

the structure from C„ to C„ , . • Let additional cartesian displacements of 

N N+l 

the material point in moving from C„ to C„ , . be Au. ; thus x’! +I = X N + Au. 

° N N+l 1 l l l 

are the new coordinates of a material point in C„ , . . The incremental Green 

N+l 

strain A g. . from C„ to , , , with reference to the metric in C„ can then be 

N N+l N 

written as , 


dAu. dAu, 


. i — i + — 1 1 + 1 r 

S ii 2 L H ,,K J 2 L 


- Siu K SAu k! _ 


" 1J "-ax'! ax'? 

J - 1 


N 

SX. Ax\ 
1 J 


N J “ Ae ij + A ^ 


* 


(3.1) 


During the motion from C N to C N+1 , the state C N is treated as one with 

"initial stresses". The new stresses in C„ , . due to the additional incre- 

N+l 

mental loading, will be represented by the symmetric second Piola-Kirchhof f 
N+l 

stress tensor , which is measured per unit area in C . 

ij (N) N 


s H+ ;„, - t h . + a* s . . 

lj(N) Lj 1J 


(3.2) 


Assuming the material has yielded, the relation between the incremental 
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^ *!* 
stress (AS..) and the incremental strain A g.. will be written as 
ij 


* t N * 

h S ij E ijkX (T ij^ h 8 k a 


(3.3) 


where is the current constitutive property, as modified by plasticity, 

N 

and is a function of current true stress, T... Such a constitutive relation, 

ij 

based on the well-known Huber -Mis es-Hencky yield criteria, Drucker's 
normality condition, and a Prager-Ziegler type kinematic hardening law has 
been derived for the present plane stress conditions, following the procedure 
described in Ref. 13. 

Based on this, it can be shown that the variational principle governing 

the equilibrium of state C.,,. is 6An = 0 where, 

N+i HD 

P 

An (Au.,Av.,T .) = X / \ J 'fr 1 ?. + ^ A*S . .) A*e.. + j t N .Aii .Au . 

HD i’ l Li/ 4^ l ^ L\ ij 2 iy xj 2 ij k, l k,j 

m 

- (F 1 ? + AF . )Au .Ida - f T.(Au. - Av.)ds - f (T. + AT.)Av.ds( (3.4) 
x x iJ ^ Lx x i' J v i l x J 


where 


domain of the m finite element. 


= boundary of the m finite element. 


= a portion of 30^ where tractions are prescribed, 

= true (Euler) stress in , measured/unit area in and in the 
fixed cartesian metric, 


F. ,T. = body force and surface traction, per unit volume and per unit 

1 ‘ area, respectively, in C^, 

* 

The S tress -increment AS., is commonly referred to as the Truesdell stress- 
increment. This stress-raie is objective [15,163. 
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(3.5) 


. *. . 

Ae ij 


, r 5 Au . 

1 1 i 

2 L -,_N 
3x. 


3 Au ... 


3x 


N J : 


AS 


ij 




* 

A e 


k£ 


Au. = arbitrarily assumed displacements, in each element, that need 

not satisfy interelement compatibility a priori. In elements 
surrounding the crack- tip, displacements corresponding to the 
Hutchinson-Rice-Rosengren singularities, for hardening materials , 
are included. 


Av^ = independently assumed displacements at the element boundary, 

3 0 , which inherently satisfy inter-element compatibility 

m 

criteria. 


T . = Lagrange Multiplier to enforce the compatibility condition, 

Au. = Av. at 30 , and 

l l m 

AF^jATL = prescribed increments of body forces and surface tractions, re- 
spectively, measured unit volume and unit area, respectively in 

V 

The Euler equations corresponding to 6Arr (6Au., 6Av., 6 t . ) = 0 lead to, 

1. 1 LI 


A*S, . . + (t^.Au. .) . + AF, + -(tJ 1 . . + F, \ = 0 

kj,j v tJ k l kj,j kJ 


( 


N * \ N 

t. . + A S, .in. + t. .Au. .n. = T_. at 30 
kj kj/ j ij k,j l Lk m 


(3.6) 

(3.7) 


T. + AT. 
x x 



at S a 

m 


(3.8) 


and 


Au . = Av . at 3 0 

li m 


(3.9) 
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In the above, a comma followed by an index such as i, is meant to designate 

N 

a partial differentiation with respect to the cartes ian coordinate , X^. 

Eq. (3.6) refers to the equilibrium of the total second Piola-Kirchhof f 
N+l 

stress S . . /X , x in . . . If state <1, was in true equilibrium, the last two- 
ij (N) N+l N 

bracketed terms in Eq. (3.6) would be equal to zero. However, because of 
the .inherent numerical errors in the incremental solution process, the 
state C N may not be truly equilibrated. Thus, retaining the last two terms 
in Eq. (3.6) leads to an "equilibrium-check" iteration process, similar to 
the one described in detail by Hofmeister, et al [17]. Eq. (3.7) states that 
the tractions derived from the assumed incremental interior displacements 
match the independently assumed boundary tractions, T . . Finally, Eq. (3.9) 

Jjl 

is the statement of inter-element displacement compatibility that is enforced 

in the present method, by means of Lagrange Multipliers, T .. 

Li L 

Because of the advantages of the above described hybrid-displacement 
model, and the convenience of the conventional finite -element displacement 
model, a combination of both the methods is used in the present formulation. 
Thus one can visualize the domain of the cracked structure to be divided 
into two regions, (a) a small region near the crack-tip where the singular, 
near field solution is predominant, and (b) a region away from the crack-tip 
where the effect of the singularity is not dominant. In the present calcu- 
lations, a hybrid-displacement model is used to derive the stiffness properties 
of the. near- tip sector elements, and the conventional displacement model is 
used to derive the stiffness properties of the far-field, eight-node 
isoparametric elements. 

Consider first the development of the properties of the near-tip, 
circular-sector shaped "singular" elements, wherein, the three Independent 
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variables Au., Av., and T t . are assumed as follows 
1 1 Li 


* 


Au. = A 5, Av. = L Aq; T . = R o 


(3.10) 


where £, and a are unknown independent parameters, and ^Aq are increments of 
nodal displacements during the generic load step, -* C^ + ^. The functions 
A and R are arbitrary. However, L are functions ar the boundary of the 
circular-sector element such that they uniquely interpolate for Av at the 
boundary in terms of the relevant Acj at the boundary, and in addition, en- 
sure displacement-compatibility with the surrounding 8-noded isoparametric 
quadrilateral elements (Fig. 2). 

For far-field regular elements, as mentioned earlier, the conventional 
compatible displacement finite element model is used. The incremental 
energy functional corresponding to this model can be expressed as. 


m 

C 


tii (Au) = , 

CD 1 -f-i.i t 

m=p+l U 


l (2 E ijkl Ae kl Ae ij + 2 T ij Au k,i Au k,j 


m 


- AF.Au. Ida - f AT.Au.ds + P (t^.Au. . - F*?Au.)da 

11/ J 11 ^ tj i,J 1 1 

o m 

m 

- f* T^Au.ds 
•l 1 r ) 


(3.11) 


m 


The finite-element approximations for a far-field regular element are con- 
sidered as. 


^Au = D Ajcj 


(3.12) 


The usual notation, (~) under a symbol, is used to denote a column vector, 
and (&) under a symbol is used to denote a matrix. 
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from which the strains and deformation gradients are derived as. 


~ = I A -s A Sk,i -gk£ 


(3.13) 


For purposes of convenience of notations we consider that, out of a total of 
M finite elements, m = 1, p are the so-called singular elements and m = p+1 
. . .M are the far-field regular elements. 

The assumptions for the field variables as in Eq. (3.10) for elements 


m = l,...p are substituted into Eq . (3.4) and the first variation of rr 
with respect to the parameters a and p is set to zero, to obtain, for 
elements m = 1, p, that 


HD 


a = 


£ - 


p" T (H +C )P -1 G Aq + p" T C 


~e 


-T 

P F. 
~4 


P _1 G Ag 


(3.14) 

(3.15) 


where , 


Ae = W P , Au, . = W, P 
~ ~ * ~k, i * 

H = P W T E WdA 
m 

P = f R T A ds 

~ s' n ** ” 

m 

G = f R T L ds 

« 5 J Q « « 

m 

C “ J* ' ( w i T ” • w i + W 9 T ” • W )dA 

~g J W «rj «1 «2 ~ij ,£.1 

m 


(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 


and 


C = f W T t N dA 




n 


m 


(3.21) 
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Substituting for a and £ from Eqs. (3.14 and 15) into Eq. (3.4), one can 

express in terms of Acj only, as, 

P 

‘ I 2 AaV + K )Aj - (3.22) 

m-1 ~ 

where , 


and 


K = (P-V H(P _ 1 G) 

K = (P _ 1 G) T C (P _ 1 G) 

= £2 = S b T h Z ds 

S a~ 

m 


'Sc = <| 1 « T <£4 


c ) 

~e 


(3.23) 

(3.24) 

(3.25) 


(3.26) 


Likewise the assumptions for the field variable in Eq. (3.12) for 
elements m = p + 1 --M are substituted into Eq. (3.11) to obtain Att^, 
for elements m = p + 1 M, as, 


A " CD * Ji 1 4 3 T <lf + £g) As - a S T £2 - As T «4 + & <3 - z7) 


where , 


K * 1 = J B T E fc B dA, 


n 


m 


Q 


k k = J (W^T N .W 1+ wj t“ w ) 

=ag ^ «1 s»ij «1 «2 j«2 


T N 


dA 


m 


c R - r b t t b dA 


T _N 


Q 


m 


(3.28) 


(3.29) 


(3.30) 
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(3.31) 


F~ = P D T AT ds 

f*2 0 ^ ~ 


m 


and 




(3.32) 


m 


Combining Eqs . (3.22) and (3.27), the total energy functional for the 
system, Arr, can be written as, 

A " = Ji [1 A a T « + £g> A a - a 4 T 'a - A a T ^ c _ 


M 

+ 2 ., 

m-p+1 


B A aV 


+ ^) A a - A a T “ r - A a lA Scj ri -' n > 


where , 


» n R T-.R J A/-\R — T»^ -L. ^R 

= £ 2 ; and ^Sc ~4 + £e' 


(3.34) 


By expressing the element incremental nodal displacements A<j in terms 

of independent generalized global displacements Aq , and using the condition 

* 

of stationarity of Att with repect to Aq , one obtains the final expression. 


l k( vv] K )“ +1 ■ -hw) + KJ 


(3.35) 


where 





r. *-)N+l 

fA i ] » 


tAQ %- p n )1 


= generalized nodal displacements in state C N (with loads P ) 

= the tangent stiffness matrix at state C„ which includes 

N 

the effects of initial stresses, displacements, and the 
influence of plastic yielding at state C^, 

= incremental displacements of the structure from C„ to C„ , , , 

N N+l 

= incremental loads, and 
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f £0 } = residual nodal forces to check the equilibrium of state C„ . 

ec N 

The solution of Eq. (3.35) with a Newton-Raphson type equilibrium correction 

iteration, is essentially similar to that presented by Hofmeister et al [17], 

and further details are omitted. Finally, we note that in order to carry 

our the incremental solution from to C^ + 2 > t * ie second Piola-Kirchhof f 

N+l 

stresses are converted to Euler (true) stresses in using the 

relation. 


N+l 
r . . 
ij 


Det 


, ax™ ax™ 

1 1 c N+1 

axfi a=4 ax" k * w 

t J 

1 


(3.36) 


4. Details of Finite Element Field Assumptions 
4.1 Circular-Sector Shaped Singular Elements 

The dominant singularities in strains and stresses at the tip of a 
crack in a monotonically loaded two dimensional body of a nonlinear material 
obeying a power-hardening law of plasticity were analysed by Hutchinson [ill 
and Rice and Rosegren [12]. (Hereafter these are referred to as the H.R.R 
singularities.) Their analyses were based on a deformation theory of plasticity 
however, to an approximation, these are considered to be valid, in the pre- 
sent work, when a flow theory of plasticity is used. Further, the analysis 
in [11,12] are valid only for stationary cracks. However, the nature of 
singularities at the tip of a moving crack in an elastic-plastic solid is 
still largely uncertain. The existing attempts in the literature to obtain 
analytical solutions to the problem of a "steadily" moving crack-tip in 
perfectly plastic materials, by Chitaley and McClintock [ 183 and Rice [19], 
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show a logarithmic strain-singularity at the crack- tip. However, one of 
these solutions [18] has been recently argued by Broberg [20] to be in error. 
Due to these reasons, and for lack of criteria to. define "steady" state 
conditions, a priori, the H-R-R singularities [11,12] are allowed, in the 
present analysis, to be present at the crack-tip which grows under fatigue 
loading. This may, however, be viewed as an approximation in the general 
context of the finite element method, in the sense that the hypothetical 
exact solution is approximated by a set of assumed basis functions. It is 
also noted that in spite of this approximation, the crack surface deformation 
profiles, after fatigue carck extension, to be discussed later, do suggest 
a possible change in the nature and order of strain singularities near the 
growing crack-tip. Finally, we note that if the nature singularities near 
the advancing crack-tip both in the "transitory" as well as "steady state" 
conditions is clarified analytically, it is, in principle, possible to 
affect the appropriate changes in the present finite element modeling. 


The H-R-R solution for the dominant solution near the crack- tip for 

' £> " 

stresses, strains, and displacements, respectively, can be expressed as: 



- 1 / +1 

0 . . k r CT . . ( 9 ) 

ij a ij T 

(4.1) 


1 n/ ^n+l- . . 

e ij~ k e r S ij (cp) 

(4.2) 

and 

1 / ., 

. n+1- , . 

u . ~ k r u. ( 9 ) 

1 e 1 T 

(4.3) 


where (r and 9 ) are polar coordinates centered at the crack-tip; n is the 
power hardening exponent in the Ramberg-Osgood relation between the equivalent 


stress a and equivalent plastic strain e , viz., e ~ a ; k and k are 

p p U 6 


- n 


constants; and o\^(cp), 6 ^( 9 ) anc * u ^C c P) are functions of the angular co- 
ordinate 9 (see Fig. 1). 
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The above dominant singularities are embedded in the sector shaped 
"singularity" elements surrounding the crack-tip. More over these singular 
functions are augmented, in the singularity elements, by a large number of 
arbitrary polynomial functions for displacements and strains. Due to this 
reason, the singularity elements in the present analysis can be fairly 
large; in the numerical examples presented later, these elements are, in 
general, chosen to be about 10 ^ times the crack length. 


4.1.1 Displacements in the Interior of Singular Element 

Based on the above reasoning, the incremental displacements in the 
singularity elements are assumed as: 


Au = P,r + 2j3 o r0 + P o r0 2 + P c r 2 + P £ r"0 + P.^cos 9 n+1 


r r l 2 


, 0 1/ n+l„ ^ a 1 n+l 0 2 , Q . 1/ n+l 

+ p^^coscpr 9 -f P 12 coscpr 9 + p^sincpr 


13" 


1 / ,, 1 / 9 

+ P^^sincpr 0 + P^^sincpr n 0 


(4.4) 


Au 0 = P 2 r + P 4 r0 + P ? r0 2 + P g r 2 + P 9 r 2 0 - P 1() sincpr n+1 


1 / , . 1 / . . « 1 / . . 

- f3^sincpr n 0 - £^ 2 sin9r n 0 + jB^coscpr n 


1 / 4-1 17 2 
+ P 14 coscpr n+i 0 + p^coscpr 9. 


(4.5) 


15 


where 0's are undetermined parameters, cp is the global angular coordinate; 
and 0 is the angle measured from the symmetric axis of the sector element, 
as shown in Fig. 1. The above are augmented by the three rigid body modes 
for the element. Thus, using the matrix notation we write. 


Au = A £ 

as in Eq. (3.10). The strain increments are derived as. 



Ae 


SAu /Sr 

> 


rr 


r 


{Ae} = - 
*■ ^ 

CD 

CD 

<3 

► ” 4 

An /r + (1/r) 

3Auq/S9 


l“*«. 


(1/r) SAu r /S0 

+ SAu Q /S - Au- / r 

t/ T “ ^ 


(4.6) 


(4.7) 


4.1.2 Boundary Displacements for "Singularity" Elements 

Along the circumference of each of the present sector shaped singular 
elements, where 3 nodes are located as in Fig. 1, each of the displacements 


Av and Av_ are assumed in the form: 
r 0 


Av. = a, , + a_ .0 + a_ .0" 
x lx 2x 3x 


(i = r ,0 ) 


(4.8) 


where a^,' a 2i anc * a 3 i are ex P resse< ^ in terms of the respective displacements 
at the 3 nodes along the circumference. Since this displacement varies 
quadratically at the circumference, continuity of displacements at the in- 
terface of the present singular elements and the surrounding 'regular' 
quadratic isoparametric quadrilateral elements is inherently satisfied. 

Along each of the two radial boundaries of the singular sector element, 
on which 4 mode (including the crack-tip as a node) are located, the dis- 
placements are assumed as, 

1 / 


Av. = b + b„ .r + b„ .r^ + b. .r n+ ^ (i = r,6) 
x lx 2x 3x 4x 


(4.9) 
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wherein the coeffecients b^...b^ are expressed in terms of the respective 
displacements at the above mentioned 4 nodes. It is seen that the above 
boundary displacement field along the radial lines is inherently compatible 
at the interfaces with the adjacent singular elements. 

Thus there is interelement displacement compatibility not only amongst 
the singularity sector elements themselves, but also with the surrounding 
regular isoparametric elements. 

4.1.3 Boundary Tractions for Singularity Element 

The boundary tractions T ^ for the singular element as defined in Eq. 
(3.10), which are mathematically the Lagrange Multipliers in the hybrid 
variational principle as in Eq. (3.4), are assumed as: 


T_ . = v . a . . 
Li l 


where a are derived from an Airy Stress function, ■if , as 


°rr = (1/r) d ' 1 ' /dr + (1/r2) S 2 ^ 02 
°q 0 = c) 2 1>/3r 2 

= (l/r^)d^/99 - (l/r)3^iJf/Srd9 


(4.10) 


(4.11) 

(4.12) 

(4.13) 


and in Eq. (4.10) denote the direction cosines of a unit outward normal 

to the boundary. The Airy Function is chosen such that the stress field 
-1/n+l 


contains a r 


type singularity, as: 


2„2 


_2„3 


4 = CK.r + 0 + a r 0 + a, r 0 + a r + a v 0 

I i 3 4 5 b 

. 3 q 2 . 3 fl 3 4 4 fl 4 fl 2 

+ 0 + a g r 9 + o^r + c^r 0 4- c^r 9 


. 4 n 3 , (2n+l/n+l) . (2n+l/n+l) . (2n+l/n+l)„2 

+ a u r 9 + a 13 r v + e + 'Q 

(4.14) 
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where a's are undetermined parameters. We note that the element boundary 
traction field as chosen above is self-equilibrated. 

4.2 Displacement Field Assumption for Regular Elements 

As noted earlier, the region of the cracked structure outside the 
above discussed singularity core region is modeled with conventional iso- 
parametric displacement elements. An eight node isoparametric quadrilateral 
finite element, as shown in Fig. 2, is used for this purpose. The trans* 
formation from cartesian coordinates of the element, with general quadratic 
curves as boundaries, to curvi-linear coordinates § and 7] , as indicated in 

Fig. 2, is achieved through: 

8 8 

x = ) N.(S,1))x. ; y = 

i=l 1 1 

where, the functions 17^(5 ,T|) are defined as: 

N.(5,T1) = (1/4) (1 +SS.X1 +T1T].)(1I. - - 1) (4.16) 

for corner nodes i = 1, 3, 5 and 7 as indicated in Fig. 2, and 

N ± (5 ,Tl) = (1/2) (1 - 5 2 )(1 +T]T1.) (4.17) 

for nodes with 5^ = 0, and, finally, 

N i (5,Tl) = (1/2) (1 - T} 2 )(1 + 55.) ’ (4.18) 

for nodes with 7)^ ~ 0. The displacement field for the regular element, 
as defined in Eq . (3.12), is expressed as: 

^ Au 

Uu) - | *; -D Aa (4.19) 

Au 

y 


A 


N ± (? ,T))y. 


(4.15) 
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where , 


[D] = " 

i r — 


2x16 


N : ; 0; N 2 ; 0 ; — N g ; 0-, 
0 ; N x ; 0 ; N 2 ... 0 ; N g J 


(4.20) 


and L Aq J = l_Aq^; Aq 2 ... Aq g _l. FromEq. (4.19), the element strains 


are 


defined as: 


Ae 


r XX , 

{Ae} = -j Ae . .!■ = B Au = B D Aq 
~ l yy J 

A\ 

xy 


(4.21) 


where [b] = 


3N . /SX 

l 


0 


0 dN^/ 3 y 

5n . / ax aN . /ay 

1 X 


(4.23) 


Finally, we note 


r aN./ax, , aN./a^ 

{ 1 } - [j ]- 1 i 1 } 


SN^/ay 


aN i /air 


(4.24) 


where [j] is the Jacobian of the transformation used in Eq. (4.15). 


5. Elastic-Plastic Constitutive Relations 

In the present analysis, based on an incremental theory of plasticity, 

* 

the piecewise linear relation between the stress and strain rates is taken 
as : 


*t 

da. . = E . . de . 

ij xjkx k x 


(5.1) 


*t 

where E.., „ is the current tensor of elastic-plastic constitutive law; and 
ljkx 

the total strain rate is assumed to be a sum of elastic and plastic parts, 


First, we develop the constitutive relation for a small-displacement theory 
of plasticity, and later generalize it for the case of finite deformations. 
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as : 


de, „ = def . + def 

kX kX kX 


(5.2) 


The elastic strain increments are assumed to be given by the Hookes' law as: 


da. . = E . .. . de . 

lj ijkX kX 


(5.3) 


where E^^x is current elasticity tensor. In the classical rate theory of 
. plasticity, which is used in the present, the elastic-plastic behaviour 
of the material is fully described by specifying: an initial yield criterion; 

a flow rate; and a hardening rule. The elastic-plastic constitutive tensor 
E i jkX. obtained by specifying each of the above, as in the following. 

(i) Initial Yield Criterion 

The well known flow theory of Huber-Mises -Hencky is used in specifying 
the initial yield criterion, which states: 


f(J 2 ) - <f> ’ 0 


(5.4) 


where a„ is the yield stress of the material and a. . is the deviatoric 
Ys 7 ij 

stress as defined by 



6 . . 



, 1 if i 

| 



6 . . 

ij 


1J ''0 if i f j 

For the special case of plane stress conditions, 
as described by Eq. (5.4) reduces to 

2 2 2 
f(a. ) = a + a + a a +3 a 
ij xx yy xx yy xy 

and for the case of plane strain conditions. 


(5.5) 


The initial yield surface 


-a 2 =0 (5.6) 

ys 
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2 2 2 2 2 
f (a . .) = a +c + a - a a - a a - a a +3a -a =0 
xj xx yy zz xx yy xx zz yy zz xy ys 

(5.7) 

(ii) Flow Rule 

The well known Drucker's normality rule, which states that the plastic 
strain increment is always orthogonal to the yield surface, or 


ds i j ‘ dX dr‘ < 5 - 8 > 

ij 

where dX is a proportionality, is used. If the plastic strain vector is 
superimposed on the stress or space (as best depicted in Fig. 3), for small 
increments of stress and strain the relation between the stress and the strain 
increments and the yield surface will be 


< da ij - alfj = 0 


(5.9) 


From Eqs. (5.8) and (5.9) the proportionality constant dX is determined as 



(5.10) 


where c is a constant. 

Substituting the above relation for dX into the Eq. 
a flow rule. 



(5.8), now we get 


(5.11) 


For the simplest case of uniaxial tension or compression Eq. (5.11) 
reduces to, 


de 


P 



or c 


dc 


de 1 


(5.12) 
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Therefore in the case of the uniaxial loading the parameter 'c' is 


the slope of stress-plastic strain curve. Or 


o o dCT 

= 4 *•' =| — ^ 

3 3 d€ P 

eq 


(5.13) 


where H y is the slope of the equivalent stress-equivalent plastic strain 

curve, and o and de p are defined, respectively, as 
eq eq 


a 

eq 


-J 


3 / / 

4 o. .a. . 
2 ij ij 


(5.14) 



j - 

J 3 


de? de P 
i-J ij 


(5.15) 


(iii) Hardening rule 

Prager's [21] kinematic hardening low, as modified by Ziegler [22] 
which appears to adequately model the Bauschinger effects that are generally 
observed in most of metals under cyclic loading is chosen. In this case, 
the subsequent yield surface is represented by. 


f (cj . . - a . . ) = 0 
ij iJ 


(5.16) 


Ziegler's [22] modified kinematic hardening rule states that the translation 
of the yield surface, da^, i- s always in the direction of the vector cr_ - 


a . . , ie . , 
IJ 


da . . = d |i (O . . 
ij ij 


a. .) 
1 J 


(5.17) 


as shown in Fig. 3. For an infinitesimal increment of loading the vector 
do . . - da., must be orthogonal to the outer normal to the surface; thus 



The scalar quantity d (J. is obtained by substituting Eq. (5.17) into (5.18) 
to get, 


d p = 


df J 
3 - — dcr. . 
da . . 11 

AJ 1_ 

df 


0) 


(5.19) 


k X kX' da. 


kX 


Therefore the rate of the yield surface translation is given by 


do . 


ij 


(cr. .-a. .) r —~- da 
rj ri da mn 
mn 

/ \ df 

(CT kX" Q 'kX ) da 


(5.20) 


kX 


Upon substitution of de^\ from Eq. (5.11) into Eqs . (5.2 and 5.3), the relation 
between the stress rate and the total strain rate can be expressed as 
follows : 


a. for plane stress case: 

2 


dcr 


xx 


da 


>= 


yy 


da 


xy 


IV 


1 Ev 

q * 2 

b 1 -v 


S 1 S 2 
S ’ 


S 1 S 3 


S 2 S 3 


1-v 


sym - 


2 ( 1+v ) S 




de 


xx 


de ^ 

yy 


2 de 

xy 

v y 


(5.21) 


where 


O = E ( df df N | 

S 1 . 2 Vda + v da ) 

1-v xx yy 


- E ( df , df s , 

. 2 V v da + d 0 ) 

1-v xx yy 


S„ = 


E df 


3 1+v da 


xy 
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s = 


CQ + s x ^ + s 2 + 2s 3 g- 

xx yy xy 


o - (^) 2 * (#-) 2 * <#-? 

XX yy xy 


_5f_ 

5a 


xx 


= 2 (a - a )-(a -a ) 
xx xx yy yy 


5f 

5a 

yy 


= 2 (a - a ) - (a -a ) 

' XTXT in?' ' TTTT w 7 


yy yy 


XX XX 


= 3 (a -a ) 
5a xy xy 

xy 


(5.22) 


<F 


b. for plane strain case : 

.2 


XX 


■% da >- 

yy 


da 


xy 


VE 


S 1 S 2 


(l-K>)(l-2v) S » (1+v) (l-2v) S 


(l-v)E _ _2 


(1+v) (l-2v) 


2 

S ’ 


S 1 S 3 


S 2 S 3 


sym - 


2 (1+v) S 


f de ^ 


xx 


J de yy ^(5.23) 


2de 

v xy y 


where 


S 1 (1+v) (l-2v) L 


xx yy 


s„ = 


2 (1+v) (l-2v) l 5a 


i^+a-v 3£ ' 


xx 


5a J 

yy 


s„ - 


E 5f 
1+v 5a 

xy 


S = ■, 


Q = 


j" are defined to be same as in the plane stress case, 
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[(ct -a ) - (c -o' )] 

oc 2 xx xx yy yy 

XX 


= J [(a - a ) - (a -a)] 

da 2 yy °yy xx xx ; 


= 3(o -a ) 
do xy xy 

xy 


(5.24) 


As noted earlier, the above relation, do. . = E . , de, . was derived for a 

xj xjkX kZ 

small displacement theory of plasticity. 

In generalizing the above relation to finite deformations, one needs to 

make a choice for an objective stress rate and an appropriate strain rate; 

and E fc „ can be used as constitutive-relation tensor relating the chosen 
ijki. 

objective stress-rate and its conjugate strain rate. In [15,163, and 
elsewhere, the advantages of using the corotational rate of Kirchhoff stress 
[which is J times the true stress T^], (which stress rate is physically 
the stress increment measured in a coordinate system that is instantansous ly 
coincident with the current carstesian system in the updated configuration 
Cjj» but is rotating with respect it at the same rate as the spin suffered 
by the material fibers in going from -* C^ + ^) , in tbs elastic-plastic 
constitutive law for finite deformations have been discussed. If AK^ . is 
the corotational rate of Kirchhoff stress, it is shown in [15,16], that it 
is related to the Truesdell stress-increment A of Eq . (3.2), through: 


* j * n N * 

A S . . = AK. . - A g ■ t" . - t‘ .. A g. . 
ij ij lk kj lk & kj 


(5.25) 


N * 

where T, . and A g, . are as defined earlier. Now, as discussed above and as 
kj 6 kj 

elaborated in [15,16], a generalization of the elastic-plastic law can be 


written as: 
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(5.26) 


J *t * 

M ij = E iju A 

By substituting Eq. (5.26) in (5.25), the tensor th e relation 

A * s ij * < 5 - 27 > 

is found, and is used in the developments discussed earlier, ie . , Eq. (3) 

on words. Finally, we note that in the above generalization for finite 

deformations, the initial as well as subsequent yield functions are expressed 

in terms of the true stress, T... 

ij 

5.1 Knee Correction 

Once a point has yielded, the state of stress represented by its location 

in the stress space stays on the yield surface until unloaded. The subsequent 

^ f 

loading condition of the point is determined, by checking-^ — dT„, as 

Tij 1J df 

loading, neutral loading, and unloading depending on whether ^ — dT„ 

T ij 1J 

is positive, zero, or negative respectively. Howevet, some complications 
arise when an elastic or unloaded part of the structure near yield becomes 
plastic with the next increment of load. A technique, herein called "knee 
correction", is used for dealing with each integration point in this transi- 
tion region, of elastic and plastic behavior, and is briefly described in 
the following. 

First an estimate of strain and stress increments which made the point 
yield is made by assuming, 


f (t . . + vAt . . - a..) = 0 
ij ij ij 


(5.28) 


where At., is the total increment of stress obtained for at the present in- 

crement of loading, and y is the proportionality constant so that the increment 

yAT . . causes yield at the point. Solving Eq . (5.28) for y, yAT . . is treated 

i J ^ J 
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as if it acts elastically, then subsequently the behavior is assumed to be 

* 

elastic-plastic, with the remaining strain increment, (l-y)A g... 

6. Finite Element Modelling of Crack Growth 

The steps in the finite element simulation of crack growth in the 
present procedure may be described as: (i) geometrical change in the crack 

surface boundary; (ii) translation of the crack- tip singularities to the 
advanced crack- tip; and (iii) release of surface tractions on the newly created 
crack surface. 

The change in the crack surface boundary is made by translating the 
whole set of crack- tip core elements, as shown in Fig. 4, by arbitrary 
distance Aa in the direction of intended crack extension; thus the new 
crack-tip node which is designated by the center of the sector-shaped core 
elements need not be coincident with any previously existing finite element 
node before extension. Thus even though the fixed boundary in the uncracked 
ligament of the structure is changed, the constraining condition of the nodes 
need not be altered. Elements immediately adjacent to the core must be re- 
adjusted to fit to the translated core. This process of translating the 
5 

core mesh also moves the embedded singularity in the elements to the new 
crack tip area, leaving no singularities but large deformations and strains 
in the wake of advanced crack-tip. All the 5x5 Gaussian data points in each 
of the translated core element (and also the 5x5 points for the conventional 
elements) may generally not coincide with those before translation, at which 
plastic history data such as current stresses, plastic strains, plastically 
dissipated work, yield surface translation, etc., are available. Therefore 
the data at points in the new mesh are estimated by linearly interpolating 
data on four Gaussian points in the old mesh that are nearest to the point 
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under question in the new mesh. The simple but cumbersome mathematical de- 
tails of this interpolation and smoothing process are omitted here for the 
sake of brevity. With the fitted plastic data and the new element geometry, 
element stiffness matrices are recalculated for the core elements as well 
as for the surrounding rearranged elements and the global stiffness is 
appropriately modified. Subsequent equilibrium check iterations using the 
new stiffness of the structure correct fitting errors, if any, in the 
plasticity data in the new mesh. At the same time, the tractions over the 
distance AB (Aa as shown in Fig. 4) are incrementally removed, with equilibrium 
check iterations being used at each step, to create a new traction- free 
crack surface of length Aa . The finite element simulation of crack exten- 
sion by the desired amount, Aa, is now completed. 

7. Analysis of Fatigue Crack Growth Under Mode I Cyclic Loading 
7.1 Description of the Problem 

Through out the series of the present elastic-plastic finite element 
analyses of fatigue crack growth under Mode I type cyclic loading, a thin 
rectangular plate with a central crack and under uniform tensile stresses, 
in a direction normal to the crack-axis, at the edges of the plate, is con- 
sidered (See Fig. 5). The dimensions of the plate are: half width w = 230mm; 

and half-crack length a Q = 27.3mm, respectively. The material is considered 

to be a 2024-T3 Aluminum alloy, whose mechanical properties are characterized 

2 2 
by: yield stress, O = 350 MN/m ; and Young's modulus, E = 70,000 MN/m . 

ys 

The material is assumed to be elastic-perfectly-plastic . It is noted that 
the above problem definition is identical to that used by Newman [8]. The 
plate is assumed to be in a state of plane stress. 
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Because of the symmetries of geometry, applied loading, and material 
homogeneity, only a quarter of the cracked plate is analysed. Fig. 5 shows 
the finite element breakdown that is used. A total of 6 sector-shaped sing- 
ularity elements near the crack- tip, and 43 conventional quadratic isoparametric 
elements are employed. Some of these isoparametric elements are 6 noded tri- 
angles, while the majority are 8 noded quadrilaterals (See Fig. 5), It is 
seen that the total number of nodes in the finite element mesh for the 
quarter-plate is 171, with a total number of degrees of freedom of 311. 

The radius of the sector-shaped singularity elements is chosen as 
p = 2.8mm; ie . , p/a Q = .103. While the crack-extension per cycle of loading, 

Aa, can be arbitrary (ie., not related to the finite element mesh size) in 
the present analysis procedure, it is chosen to be Aa = 0.14mm in the present 
series of computations. 

7.2 Techniques to Minimize Computational (CPU) Time 

Firstly, we note that the near-tip elements used presently are of the 
order 10 ^ times the semi-crack length; and the total number of algebraic 
equations for the above problem are only 311. 

In the present procedure, a tangent modulus (stiffness) approach is 
used in each increment and in each iteration in the respective increment. 

Thus a faster convergence is obtained in the process of iterations of 
equilibrium correction, etc. However, it is noted that only the stiffness 
matrices of the plastic portion of the structure need to be changed in the 
present tangent modulus approach, whereas, those of the elastic portion re- 
main fixed. 

For the presently considered levels of applied far-field tension on the 
specimen, only "small-scale" yielding conditions prevail near the crack-tip. 
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A typical plastic-zone size near the crack-tip for the presently considered 
class of Mode I problems is shown in Fig. 6a, being superimposed on the finite 
element mesh. It is seen that the plastic portion of the structure (designated 
"Region P") is considerably smaller than the elastic portion (designated 
"Region E") . Thus only the stiffness matrix of region P need to be changed 
in each load-step and each iteration in the present procedure. 

However, the total number of times, say N, that the combined stiffness 
matrix (for Regions P + E) must be inverted in the course of analyzing a 
typical problem, say the case a growing crack under constant amplitude 
cyclic loading, is S = [(no. of iterations/cycle) X (number of load increments 
per cycle) X no. of load cycles ]. For a typical problem, say 8 cycles of 
constant amplitude loading, a typical value for N can beN=4x28x8= 896, 
(ie., 4 iterations per cycle, with 28 load increments/cycle, etc.). This 
is a rather enormous amount of computing; and hence a more economical way 
of solving the stiffness equations is mandatory. 

Since the plastic-zone is a small-size, by an appropriate node numbering 
scheme, the stiffness matrix of the plastic zone can be arranged, as in Fig. 
6.b, so that it is a small sub-set of the global stiffness matrix of the 
structure (eventhough the plastic zone size keeps changing with load, an 
approximate preliminary analysis can be made to determine its size at 
maximum load). Then we use a static condensation procedure to first 
eliminate the equations corresponding to the nodes in the elastic portion, 
in the very first load -increment . Thereafter only the equations for the 
nodes in Region P need to be operated upon, in all subsequent load steps 
and iterations. Thus, in the example cited above, (N = 896), in all but 
one of the 896 solutions, the number of equations being solved is rather 
very small, and correspond to the total number of unconstrained nodal 
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displacements in Region P. This enables the present computations to be 
economically feasible. 

Also, the computer program is so arranged that the data obtained from 
computation upto the end of a given spectrum of loading can be used as in- 
put data at the beginning of a different spectrum of loading. For instance, 
at the end of a constant Hi-amplitude cyclic loading, the data is stored on 
a direct access permanent disc file and used as initial conditions for a 
low amplitude cyclic load specturm; in this process not only the case of 
Hi-amplitude loading but also the case of Hi-to-Lo block loading is solved. 


7.3 Monitoring of Crack-Closure and -Opening in the Finite Element Model 

Let us assume, that at a given instant of time (at a given point in 
the loading history), the location of the crack-tip, the location of the 
A nodes on the radial line (which coincides with the crack surface) of the 
singular sector element, the locations of all other nodes on the crack axis, 
as well the current (deformed) profile of the crack surface, are all known. 

We denote the current nodes on the crack surface as "Updated Lagrangean 
Nodes". Let us now assume that the crack-tip is now further extended by an 
arbitrary amount (Aa) and the structure is then subjected to further loading. 
We first note that, in the present development, the boundary displacement 
(in the direction of the applied normal stress) along a radial line of a 
"singular" sector element is of the form. 



1 / 


a r 


n+1 , 2 

+ br 


+ cr + d 


(7.1) 


where a, b, c, and d are expressed in terms of the normal displacement, 

Av , of each of the 4 nodes on the radial line of the sector element in its 

y 

current location. Using the above equation, and knowing, a priori, the 
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radial coordinates, (ie., r, as measured from the current crack-tip), of 
the respective nodes on the radial line of the sector element in its 
immediately previous location, one can compute the values of Av^ at the 
above mentioned "Updated Lagrangean Nodes". By adding (or substracting, 
as the case may be,) these incremental displacements to the previously 
known values, an accurate prediction of the current crack surface deformation 
profile is made. 

During the unloading part of any cycle, of the present cyclic loading 
case, at the instant the displacement (in the direction of applied tension) 
at one or more nodes on the crack surface becomes negative, further un- 
loading stopped. The computational procedure is then switched to a dis- 
placement control type, and the above negative displacements are preeisely 
enforced to be zero; thus finding the precise stress level at which the 
closure constraint on the respective node must be enforced. 

After the crack-closure is detected, the respective node(s) are constrained 
thereafter, until the restraining force(s) at the node(s) just becomes zero 
and begins to be tensile in nature. The corresponding applied stress level 
defines the crack-opening stress. 

The above processes are illustrated for the case of constant amplitude 
cyclic loading in Fig. 7 ; in which, the vertical displacement (when the 
crack is open) and the restraining force (when the crack closes) at a node, 
at which the crack closes and opens, are plotted against the applied stress 
level. It can be seen from Fig. 7 that the crack-opening stress (or 

the closure stress (7 can be most accurately determined by extrapolating 
to zero the respective nodal displacement after opening (or just before 
closure, to find a c $) > as well as by extrapolating to zero the restraining 
force at the respective node just before opening (or just after closure to 
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find a „). It can be seen that the values a and a . obtained by the 
cJL op ci 

above two sets of extrapolation procedures agree excellently. This indicates 
that the present numerical procedure of determining ct and o by monitor- 
ing crack-surface displacements and the nodal restraining forces, is very 
accurate . 

Finally, some comments on the presently observed patterns of crack- 
closure are given, before proceeding to a discussion of specific cases. In 
general, closure was noticed to occur at the node closest to. the current 
crack-tip, as indicated in the sequence of unloading steps in Figs. 8a-c. 
However, if the current crack-sruface profile is irregular, as in the case 
of Hi-to-Lo block loading to be discussed later, crack-closure may first 
occur at the node closest to the crack- tip; however, in the subsequent 
unloading step, closure may occur at a node far-removed from the crack-tip, 
as indicated in Fig. 8d. From the results to be discussed later, this 
pattern of crack-closure appears to contribute significantly to growth 
retardation and delay effects. 

7.4 Criterion For Crack-Extension Stress Level 

In the present work, a study is made to arrive at a criterion for the 

stress level, a , at which fatigue crack growth occurs. In prior literature, 
ex 

this crack-extension stress level was chosen arbitrarily. For instance, 
in [8] the crack is extended at the maximum applied stress in each cycle 
even in a general spectrum (for instance, high-to-low, low-to-high, etc.,) 
loading, where as in [10] the crack was extended at the applied stress level 
at which the restraining nodal force at the new crack- tip becomes zero. In 
the present study, for instance in a constant-amplitude (zero-to- tension) 
cyclic loading, it was found that the crack opening and closure stresses. 
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o and a „ respectively, were sensitive to the chosen a 
op c* ex 

In the present work, a criterion, a = a + p(a -a ), where 

ex op max op 

p is a constant of proportionality, is postulated; and p is obtained by 

calibration such that the calculated a correlated with that observed in 

op 

experimental studies such as in [l,2,3]. However, it is noted that this 
constant of proportionality p may, to some extent, be dependent on the 
numerical methodology employed in fatigue crack modeling itself. Thus 
the above described calibration may be considered as valid only in the 
context of the particualr methodology employed in the present work. 

Three different test cases, each with a different magnitude of constant 
amplitude (zero to tension) cyclic loading, were studied with different 
values being chosen in each case for the above mentioned constant of 
proportionality, p. The idea was to select a ' p * that yields results, in 
each case, for (a /c^^) that are in best agreement with the experimental 
results [l,2,3l for 2024-T3 Aluminum alloy, which is the material simulated 
in analysis. 

The results, for instance, for the case (o' /a = .40) and 

max ys 

(R = a . /tf - 0) are summarized as follows: 
mm max 


p 

Levelled-off cr 

op 

[a /a ] at steady state 

op max 

1.0 

115 MP 

.82 


a 


.85 

94 

.67 

.62 

79 " 

.56 

.40 

58 

.41 
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Similar results were obtained for the cases, (0 fa ) = .229 and 

max ys 

(o /a ) = .314. From these three sets of results, it was observed that 
max ys 

p~0.62 yields results for (o /cr ) that are in best agreement with 

op max 

experimental observations, which indicate that [a /a 1 at steady state 

op max 3 

is about 0.56. 

Thus when p is chosen to be 0.62, the results obtained for the steady 

state vaule of (a /a ) for three-different cases of (O /cr ) are 
op max' v max ys 

summarized as follows: 


Case No. 

cr /a 

max ys 

Levelled off a 

op 

(cr /a ), steady state 

op max ’ ' 

1 

0.229 

44 MP 

0.55 



a 


2 

0.314 

57.5 MP 

a 

0.52 

3 

0.40 

78 . 5 MP 

a 

0.56 


It is hypothesized that the above constant p = 0.62 may be used throughout 

the rest of the analysis, ie . , for cases of general spectrum loading. We 

also note that when the load level cf during any cycle is first determined, 

the number (and size) of load steps between this o and o , in the 

op max 

respective cycle, is so adjusted that the pre-chosen level of = c + 

p(a - a )J coincides with one of the load- increments in the cycle, 
max op 

We now discuss the results of analysis of each of the loading cases. 


7.5 Constant Amplitude Zero-to-Tension Cyclic Loading 


(i) The results for a , for the case of (a /cr ) = 0.4 and R = 

op max ys 

(cr . /a ) = 0, are shown in Fig. 9, for 8cycles of loading. It is ob- 
v min max & 3 

til 

served that cr reaches a "steady-state" value of 0.56 o after the 4 
op max 
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tlx 

or 5 1 cycle. It is also noted that this value for o /o (= 0.56) iq 

op max 

in reasonable accord with experimental results [1,2,3] for the same material, 
a 2024-T3 Aluminum alloy. 

(ii) From Fig. 7, wherein the opening displacement and restraining force 
at a node, at which crack-closure occurs, are plotted against the applied 
stress, a nonlinearity in the variation of the opening displacment, as the 
load approaches ff (for instance, between points D-E in Fig. 7 ) is observed; 
this can be seen to be due to plastic deformation. In the subsequent un- 
loading (from point E in Fig. 7 ) the opening displacement decreases linearly, 
thus showing the effects of elastic unloading. When crack closure occurs 
(point F in Fig. 7), a nonlinearity in the variation of the nodal-restraining 
force (between F-G) can be noticed until the applied stress in the cycle 
becomes zero; this reflects plastic compression. 

(iii) Knowing a in each cycle, we define the effective stress- intensity 

o op 

factor as : 


AK 


ef f 


= C, 


Vn(a o + NAa) 


( CT . 


max 


- a ) 

op 


(7.2) 


where C, is the finite-size correction factor for the present crack geometry 

(which was found to be = 1.017 from a finite element linear analysis of 

the crack with a = a ; and thereafter assumed to be constant); N is the 

o 

number of cycles and Aa is the crack growth per cycle. For convenience, 

the variation of AK is plotted against a, the current crack length, in 

err 

Fig. 10. It is seen that AK^^ levels-off after a few cycles to a steady 
state value. 

(iv) Fig. 11 shows the crack surface deformation profiles, for in- 
stance, at various stages of unloading during the 8^ cycle of the present 
constant-amplitude (R = o) cyclic loading. The large blunting at the initial 
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crack-tip location (a = a Q ) is observed to remain permanently. The surface 

til 

of the extended crack is observed to be fairly smooth. During the 8 
cycle, it is observed that precise crack-closure occurs only over the area 
Aa (ie., only at the previous location of the crack- tip node) even upon 
full unloading; however, it is also noted that at this stage, the opening 
of the crack faces between the points a Q and the precisely closed node is 
very small (See Fig. 11) . 


Lo-To-High Block Loading 

A two level block loading, from low to high, with C7 rnax in the higher 

level being 1.273 times the 0 in the lower level, is considered. The 

max 

maximum stress in the lower level is taken such that (0 ) T 4 (a ) = 

max Low ys 

.314. As mentioned earlier, the data at the end of 4 cycles of low level 
block loading (See Fig. 12) is recovered from a constant-amplitude test 
case, with the corresponding stress level. The following results were 
obtained : 


(i) The variation of crack-opening stress, 0 , as the cyclic loading 

progresses, is shown in Fig. (12). It can be seen that immediately after 
the step up in the level of applied stress, ff decreases by about 33% 

of its steady state value corresponding to the lower level of block loading. 
Subsequent to this, 0 increases monotonically to a steady state value 
corresponding to the higher level of block loading, within about 5 cycles. 
Prior to this stabilization, AK^^ in the higher level of block loading re- 
mains considerably higher than the steady state value corresponding to this 
stress level; thus indicating growth acceleration following the load step-up 

(ii) The variation of A K CJ . with N, with AK being defined as earlier 

ett err 

is shown in Fig. (13). It is seen that there is a jump in AK^^ immediately 
following the step-up in the level of applied loading; thereafter AK^^ 
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decreases monotonically down to its base line value corresponding to the 
higher level block loading, but doesn't go below this base line value. This 
appears to indicate that the growth rate of the crack, initially growing 
at a steady rate corresponding to the lower level blocking loading, suddenly 
accelerates and then smoothly shifts down to a steady state rate correspond- 
ing to the higher level block loading. 

(iii) The representative crack surface profiles for instance, at 
various stages of unloading at the end of the 8^ cycle (high stress) of 
the current low-to-high level block loading, are shown in Fig. (14). 

From this Figure, it can be seen that the step-up in the level of loading 
causes a blunting of the crack- tip (ie., at the location x/a^ = 1.02 in 
Fig. (14), when the step-up in loading occurs in the present finite element 
simulation). Even during the unloading at the end of the present two level 
block loading, as seen from Fig. (14), the crack-closure occurs only over 
the area Aa .(ie., only at the previous location of the crack-tip node). 


High-To-Low Block Loading * 

After 8 consequitive cycles of a high level block loading (the data 
at which point is recovered from the corresponding constant amplitude test 
case), the o is reduced by 21.4% and 8 more cycles of this reduced level 
block loading were considered. The magnitude of the applied stress in the 


high-level block was such that (cr ), . , 4 a = 0.40. The following 

6 max high ys . 

results were obtained. 


(i) The variation of the crack-opening stress O' , as the loading 
progresses, is shown in Fig. (15). It is seen that immediately after the 
step-down in the load level, no abrupt decrease in o , as was the case in 
Low- to-High loading, occurs in the present High-to-Low block loading case. 
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within 


After the load-level step down, a stayed at about 0.70 (a ), 

op max low 

the number of cycles of low-level load considered presently. It may be 
possible that, as further number of low-level load cycles are considered and 
the crack- tip grows further and eventually surpasses the plastic zone 
created during the high level block loading, the c decreases to a base 
line value corresponding to the lower level block loading. However, limita- 
tions of computer funds precluded the possibility of considering a larger 
number of load cycles at the low level. 

(ii) The AK versus N curve (plotted for convenience as AK __ vs a/a 

ett eff o 

for the present case is shown in Fig. (16). It is seen that after the load- 
level step down, AK^ff remains remarkably lower than its base line value 
corresponding to the low-level block loading; thus indicating the presence 
of a considerable retardation of growth, but no delay. 

(iii) Fig. (17) shows the crack surface profiles during various 
stages of unloading in one of the low-level cycles of the present Hi-to-Lo 
block loading. It is seen that at the stage of unloading indicated by point 
'B 1 in Fig. (17), the crack closes only at the previous crack-tip (closure 
area = Aa) . Further unloading, represented by point C, causes another node 
away from the current crack-tip to close, as seen in Fig. (17). The area 

of crack-closure thus increases as the unloading progresses. 

To understand the effects of the features of crack-closure as in the 
present case, the problem was reanalysed with the constraint of closure being 
removed on the node (as discussed above) far away from the crack- tip, but 
leaving the closure-constraint on the node closest to the crack-tip. The 
corresponding changes in a are indicated by a broken line in Fig. (15) ; 
while the change in AK^^ is also indicated by a broken line in Fig. (16). 
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These results indicate the influence of properly imposing the closure- 
constraints on nodes even far away from the crack-tip; when the considered 
loading, as the present Hi-to-Low case, causes such a type of crack-closure. 

Single Over Load 

The case of a single overload after 4 cycles of a constant amplitude 
block loading, followed by further cycles of constant amplitude (equal in 
magnitude to that before overload) was considered. The following three 
cases were considered. 


Case 

Base cr 

max 

cr overload — 
°max, base. 

a 

op - cr 

max op. base 

cr - o 

max base op .base 

1 

110 MP 

1.273 

0.151 


a 



2 

110 MP 

1.455 

0.283 


a 



3 

80 MP 

2.0 

0.681 


a 




In the above, cr , is the maximum applied stress prior- to or after 
max. base 

overload; a , , is the overload stress; o is the maximum calculated 

overload op. max 

value for crack-opening stress after overload; and o is the base line 

& op. base 

opening stress for an otherwise constant-amplitude cyclic load at level 

o , , and all these stresses are illustrated in Fig. (18a). The ob- 

max.base 

tained restuls are discussed below: 

(i) The variations of c during the load cycling, for the three different 
ratios of stress-overload, are indicated in Figs. (18a, b,c) respectively. 


In all the three overload cases, an abrupt decrease in cr (which relative 

op 
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decrease becomes more predominant as the overload stress-ratio increases) 

is noticed immediately after the single overload application. After this, 

in all the three cases, o increases again to reach a peak value c 

op ° r op .max 

before levelling off a steady-state value. This relative values of a 

opmax 

increases as the overload stress ratio increases. Also as the overload 


stress-ratio increases, the later is the occurance of this o . For 

op. max 

instance, for the overload stress-ratio of 1.273, o occurs in the 4^ 

opmax 

cycle after overload (Fig. 18a); for overload ratio of 1.455 this occurs 
til 

in the 8 cycle after overload (Fig. 18b); while for the case of overload- 
rato 2.0, a is still increasing (Fig. 18c). This implies that the higher 
the overload ratio is, the more remarkable both the retardation and delay 
effects are. 

(ii) The curve depicting the variation of the ratio (o - a )/ 

7 r opmax op. base 

(a , - a ) with the overload stress -ratio, which is drawn using the 

maxbase opbase 

above discussed 3 data points, is shown in Fig. (19). By extrapolation, 

the threshold value of the overload ratio, at which retardation effects come 

into play, is seen to be about 1.10. In contrast, Bernard et.al [ 23 report 

a threshold overload ratio of 1.3 ~ 1.4 based on a series of experiments on 

the material Ducol W30B whose .yield strength is 366 MP (, comparable to the 

a 

presently considered c = 350MP a ). It is noted however, that the present 
analysis is based on a plane-stress assumption, while Bernard et.al [23] 
note the dependence of the experimentally determined threshold value on the 
specimen thickness. 

<4 

(iii) The variation of AK with N (or a/a ) for the single-overload, 

err 0 

case, with a ratio of 2, is shown in Fig. (20), with a similar variation 
being noted for the other two overload cases considered. It is seen that 


AK experiences a sudden jump immediately after the overload, and then 
ef f 
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decreases below its base line level corresponding to a constant-amplitude 

cycling at a level; thus indicating the presence of retardation 

max 5 DasG 

and delay effects (these terms are used here in the same sense as defined by 
Barnard, et.al [23]) in this single overload case. It is also noted that the 
quantitative effects of retardation and delay depend on the overload ratio. 

(iv) The crack- line deformation profiles at various stages of unloading 
at the end of the considered number of cycles are shown in Fig. (21) for the 
case of overload ratio of 2 , while similar results were noted for the other 
two overload-ratio cases also. It is seen that the application of the single 
overload to the specimen (at the instant when a/a Q = 1.02 in Fig. 21) 

causes a large (plastic) blunting which is retained in the crack-surface 

/ 

profile even as the crack advances in further cyclic loading. When the 
specimen is fully unloaded, at the end of the cycle illustrated in Fig. (21), 
almost the whole surface area ahead of the previously mentioned location of 
blunting is noticed to close, while the crack surface area behind this 
blunting location is seen never to close. 

8 . Analysis of a Center-Cracked Specimen under Pure Mode II Cyclic Loading 

A center cracked square plate under a constant amplitude cyclic loading 

of pure shear, which is uniformly distributed at the edges of the plate, is 

analysed. Plane stress conditions are assumed. The material is considered 

to be 2024-T3 Aluminum alloy, (same as in the Pure Mode I case discussed 

earlier). The dimensions of the plate are: L = W = 140mm; a Q = 40mm. The 

maximum amplitude of the uniformly distributed shear is taken to be t = 

max 

80 MP (t /a = .23) . 
a max ys 

In the present problem, the geometry of the plate with the crack is 
symmetrical about both the x and y axis (See Fig. 22), and the external 
loading in anti-symmetric with respect to both x and y axis. 


42 



As earlier, the present material is modeled as an elastic-perfect- 
plastic material. We note also that the presently considered material has 
the same properties in tension as in compression. 

Thus, the displacement field has the following antisymmetric properties 

u x (x, y ) = -u x (x,-y) = u x (-x,y) = -u x (-x,-y) (8.1) 

u y (x,y) = -u u (-x,y) = u y (x,-y) = -u y (-x,-y) (8.2) 

where u x is the displacement in x direction, etc. Further, it is noted that 

these displacements may be discontinuous at the crack surface, -a < x < a . 

o o 

Thus, in the finite element modeling, only a quarter of the plate is modeled 
(as shown in Fig. 22) with the displacement boundary conditions: U x = 0 

along y = 0, in the uncracked ligament only; and u^ = 0 at nodes along x = 0 
The linear elastic results, based on the first load increment, from the 
present finite element analysis indicate: 

K = 0.075; K I] . = 3.777 (8.3) 

which compare favorably with the following results (obtained by using the 
finite-size correction factors of) Bowie and Neal [24] for an identical 
problem: 

Kj = 0.0 K = 3.899 (8.4) 

The fact that t 0 in the present finite element analysis is the result 
of inherent numerical errors, such as round-off and truncation, in the 
finite element analysis. 

Also, to check the numerical accuracy of the present finite element 
modeling, a second run was made wherein a half of the plate was modeled. 
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as shown in Fig. (23), instead of only a quarter of the plate. In this 


model, u^ was set to zero at nodes along x = 0, along with constraints to 

suppress the other two rigid body motions of the plate as a whole. Both 

the sets of results (ie., those when a quarter and a half of the plate, 

respectively, are modeled) correslated excellently. 

In Fig. (24a) the results for the displacements u at the upper and 

u SL 

lower crack surfaces, u and u respectively, are plotted for the linear 

u Si 

elastic case. These numerical results for u and u are identical in the 

y y 

u Si 

linear elastic case. This equality of u^ and u^ is noticed as the loading 
continues in the first cycle (crack being stationary) and when the plastic 
zone size is significant at T = 70 MPa (see Fig. (25) for the shape of the 
plasticity zone at r = 70 MPa). 

For lack of any other criteria, the crack was extended, in the present 
procedure, at 70 MPa (the maximum applied stress being 80 MPa), in the 
direction of the initial crack ie., in the x-direction. It is seen from 
Fig. (24a) that significant increase in u^ is brought about by the process 
of crack extension and hence the attendant plastic unloading; however, again, 
u^ and u^ are almost identical (to the 4 t * 1 significant digit). Thus it is 
seen that through the all stages of loading, crack extension and plastic- 
unloading, and further loading (to 8 MPa in this case) after crack extension 
the upper and lower crack faces experience identical displacements in the 
y direction, ie . , perpendicular to the initial crack axis. 

On the other hand, the displacements u at the upper and lower surfaces 
u l 

of the crack, u^ and u^ are plotted in Fig. (24b) for the cases of loading 
when linear-elastic conditions prevail (T = 31.1 MPa); when appreciable 
plasticity develops at the crack- tip (t = 70 MPa), when the crack is ex- 
tended (and hence there is plastic loading) at T = 70 MPa, and when the load 



is further increased to 80 MPa after crack extension. It is seen that the 

magnitudes of u and u are nearly identical (to the 4 significant digit), 

but with opposite sign, in all the above cases. However it is interesting 

to observe that the change (as a ratio of the respective value prior to 

crack extension, at the same load) brought about by the process of crack 

extension (and hence plastic unloading) in is much more pronounced than 

in u . 
x 


It is interesting to note that in the linear elastic case the crack- 

u H 

surface displacements u^ and u^ (see Fig 24a) exhibit almost a linear 

Oi 

variation from the crack-tip; thus indicating a lack of any r (in particular, 
Jr type for linear elasticity) component in u^ for the linear elastic case. 


On the other hand, for the linear elastic case, the tangential displacements 

u U and u^ (See Fig. 24b) do exhibit a Jr behavior near the crack-tip, 
xx 

for the present Mode II problem. Also, it is seen from Fig. (24b) that. 


u SL 

as plasticity develops, the tangential displacements u^ and u^ exhibit a 
r a (a < 1/2; and specifically ot — 1/n+l as corresponding to Eq. (4.3)) 
variation near the crack-tip. However, for the pure Mode II case, even 
in the presence of plasticity, the analyses of Hutchinson [ 11 3 and Rice 
and Rosengren [12], indicate that there may not be a r^ (a < 1/2) type 
variation in u near the crack-tip. But the present results for u in the 

y y 

presence of plasticity, Fig. 24a, are seen to contain such an r (o' < 1/2) 
type variation near the crack-tip. However, it should be noted that the 
angular variation of the singularity functions, u^(^), as embedded in the 
present sector elements (See Eq. (4.3)) are being approximated as quadratic 
polynomials in each sector element (See Eqs. 4.4 and 4.5). The fact that 

Q/ 

r (O' < 1/2) type variations in u are numerically obtained along the radial 
line of the sector element lying on the crack sruface, as in Fig. 24a, 
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suggests that the above angular variations are not being solved highly 
accurately in the present numerical method. However, it appears that these 
numerical errors are identical at i = + rr as well as at ^ = -tt , so that u U = 

y 

JL 

as in Fig. (24a). 

Finally the crack surface displacements (u U , u'S and (u U and u^) at the 

y y x x 

end of the first cycle of loading (ie., when the applied stress is brought 

back to zero) are also indicated in Figs. (24a) and (24b) respectively. 

u JL 

Once again it is seen that even after complete unloading, u and u are 

u JL 

identical in magnitude and direction, where as u^ and u^ are identical in 
magnitude but opposite in direction. 

Thus for the present material, with identical properties in tension 
as in compression, it is seen that in all cases of pure-shear type external 
loading, the upper and lower surfaces of the crack move together in the 
direction perpendicular to the initial crack-axis, whereas they slide 
past one another in the direction of the crack-axis. Thus, it appears for 
these types of materials the phenomenon of crack-closure, as observed 
experimentaly and as analysed presently in Mode I type loading conditions, 
does not occur in pure Mode II type cyclic loading. 

.However the present experience indicates that crack-closure may occur 
even in pure Mode II cyclic loading if the material has different properties 
in uniaxial tension and compression. Consideration of such materials is not 
pursued in the present work. 

Finally, the computed shapes and sizes of the plastic zone near the 
crack-tip at various stages of pure shear loading are indicated in Fig. (25). 
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9. SUMMARY AND CONCLUSIONS 


Based on the results discussed so far, the following observations and 
conclusions appear to be possible: 

(i) If one assumes that the crack- growth-rate , da/dN, in fatigue- 
loading of Mode I type, is related to the effective stress-intensity factor 
in the form: 


da / A Tr \ t 

dN " c(AK eff ) 


(9.1) 


where AK = c, V tt a (a -a ) and n is an exponent, the present results 
for 4 illustrative types of block loading indicate that: (a) growth retar- 

dation occurs in Hi-to-Low and single overload cases, and acceleration in 
Lo-to-Hi block loading, (b) significant delay effects prior to retardation 
are observed in the case a single- over load in an other wise constant ampli- 
tude cyclic loading. 

(ii) The crack-surface deformation profiels, in the different Mode I 
loading cases, viz., constant amplitude, Lo-to-Hi, Hi-to-Lo, and single 
overload, have distinctly different characterstic shapes as sketched below: 



CONSTANT LO-TO-HIGH HIGH-TO-LOW SINGLE OVER-LOAD 

AMPLITUDE 


Schematic Representation of Typical Crack-Surface Profiles in Different 
Cases of Applied Cyclic Loading. 
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Each of the above crack surface deformation profiles leads to a 
different pattern of crack-closure when the specimen is unloaded. In the 
cases of Lo-to-Hi loading and single overload, a considerable amount of 
crack-surface blunting occurs at the instant when the applied load is 
stepped up. This blunting precludes the possibility of crack-closure 
behind the location of this blunting in the subsequent load cycling (in the 
cases of Lo-to-Hi and single overload). Thus in these two cases of loading, 
immediately after the step-up in load level takes place, crack closure 
occurs only in the small area adjacent to the current crack-tip. Consequently, 
the total compressive force acting on the closed surfaces in these two 
cases, is also small immediately after the load step-up as shown in Figs. 

(26) and (27) respectively. In Fig. (26), the maximum opening displacement, 

in each cycle (ie., at maximum applied load indicated, for instance, by 

6 

point A in Fig. 26) at the node closest to the current crack-tip, as well 

as the total compressive force over the closed area of the crack upon total 

unloading in each cycle (such as at point B in Fig. 26) are shown for the 

case of Lo-to-High loading; while similar results are shown in Fig. (27) 

for the case of single-over load . In both the cases, it is seen that the 

total compressive force over the closed crack-surface- becomes a minimum 

right after the step up in the applied load; thus resulting in a lower 

value of a at which the crack-opens in the tensile portion of the sub- 
op 

sequent loading cycle. These phenomona, brought about essentially by a 
smaller area of crack-closure (due to crack-blunting that occurs at the 
instant of load step-up), may be responsible for the 'delay' effects - 
such as, the delayed retardation in the single overload case, and the de- 
layed transition of opening stress levels from the base-line value for 
lower amplitude block loading to the higher base-line value for the higher 
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amplitude block loading (see Fig. 13) in Lo-to-High loading case. 

(iii) As mentioned, growth retardation occurs in the cases of single 
overload and Hi-to-Low block loading. Observing the crack-surface profiles 
for the forementioned two cases, at the instant of complete unloading in a 
cycle, several cycles after the initial overload, as shown in Figs. (21) 
and (17), respectively, it is seen that the area of crack closure for both 
the cases is much larger in comparison with the other two cases (see Fig. 

11 and 14). Consequently, the total compressive force in the closed crack- 
surface is relatively larger in the single-overload and Hi-to-Lo cases than 
in the other two cases. Thus, in the single overload and Hi-to-Lo cases, 
after the change in applied level, the crack-opening stresses level off at 
a much higher value than the baseline value corresponding to the level of 
applied loading after the formentioned change in loading level. 

From these observations, it appears that the pattern of crack-closure 
has a significant effect on growth retardation in general spectrum loading. . 

(iv) The crack surface profiles for the Hi-to-Lo loading case, several 
cycles after the load-step down as shown in Fig. (17), possess certain novel 
features. Crack closure is seen to be possible not only at the node closest 
to the crack-tip, but also nodes for removed from it, with no closure 
being observed at nodes in between. In this case also, growth retardation 
was found to be possible. 

(v) Thus, a correlation between the process of a significant crack- 
blunting at the instant of transition in the amplitude of applied load- 
transition, and the possibility of growth retardation is found to exist, 
as summarized below. 
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Loading Case 

Blunting of 
Crack at Load- 
Transition 

Dip in a 

op 

in the Subsequent 
Cycle 

Delay 

Area of 
Crack Closure 

Retardation/ 

Acceleration 

Constant 

NO 

NO 

NO 

Small 

NONE 

Amplitude 






Lo-to-Hi 

YES 

YES 

YES 

Small 

Acceleration 

Hi-to-Lo 

NO 

NO 

NO 

Large 

Retardation 

S ingle 

YES 

YES 

YES 

Large 

Retardation 

Overload 







, The phenomenon of crack-closure was not observed in the present numerical 
modeling of a thin-center-cracked sheet (of an elastic-perfect-plastic 
material with identical properties in uniaxial tension as in compression) 
subject to external pure shear (Mode II) cyclic loading of constant amplitude. 
The implication of this in the more general problem of fatigue crack growth 
under mixed-mode cyclic loading deserves further study. 

j • 


0 


50 




























REFERENCES 


1. Elber, W., "Fatigue Crack Closure Under Cyclic Tension," Engineering 
Fracture Mechanics, Vol. 2, No. 1, July 1970, pp. 37-45. 

2. Elber, W. , "The Significance of Fatigue Crack Closure," Damage Tolerance 
in Aircraft Structures, ASTM STP 486, 1971, pp. 230-242. 

3. Rueping, J. E., Hillberry, B. M. , Mettler, S. C., and Stevenson, W. H., 
"Comparative Study of Fatigue Crack Closure," AIAA Journal , Vol. 16, 

No. 4, Apr. 1978, pp. 413-415. 

4. Fatigue Crack Growth Under Spectrum Loads , ASTM STP 595, American Society 
for Testing and Materials, 1975, pp. 281-302. (All the papers in it). 

5. Budiansky, B. and Hutchinson, J. W. , "Analysis of Closure in Fatigue 
Crack Growth". Journal of Applied Mechanics Vol. 45. No. 2, June 1978, 
pp. 267-276. 

6. Newman, J. C., Jr. "Finite-Element Analysis of Fatigue Crack Propagation- 
Including the Effects of Crack Closure." Ph.D. Thesis, Virginia Polytech- 
nic Institute and State University, Blacksburg, Va. , May 1974. 

7. Newman, J. C. Jr., and Armen, Harry, Jr., "Elastic-Plastic Analysis of a 
Propagating Crack Under Cyclic Loading," AIAA paper No. 74-366, presented 
at the AIAA/ASME/SAE 15th Structures, Structural Dynamics, and Materials 
Conference, Las Vegas, Nev. , 17-19 April 1974. 

8. Newman, J. C., Jr., "A Finite-Element Analysis of Fatigue Crack Closure", 
Mechanics of Crack Growth . ASTM STP 590, American Society for Testing 
and Materials, 1976, pp. 281-301. 

9. Ohji, K. , Ogura, K., and Okubo, Y., "Cyclic Analysis of a Propagating Crack 
and its Correlation with Fatigue Crack Growth," Engineering Fracture 
Mechanics . 1975, Vol. 7, pp. 457-464. 

10. Ogura, K. and Obji, K., "FEM Analysis of Crack Closure and Delay Effect in 
Fatigue Crack Growth Under Variable Amplitude Loading." Engineering Frac - 
ture Mechanics , 1977, Vol. 9, pp. 471-480. 

11. Hutchinson, J. W. , "Singular Behaviour at the End of a Tensile Crack in a 
Hardening Material". Journal of Mechanics and Physics of Solids , 1968, 

Vol. 16, pp. 13-31. 

12. Rice, J. R. and Rosengren, G. F., "Plane Strain Deformation Near a Crack 
Tip in a Power Low Hardening Material," Journal of Mechanics and Physics 
of Solids, Vol. 16, 1968, pp. 1-12. 


51 



13. Atluri, S. N. and Nakagaki, M. , "J-Integral Estimates for Strain-Hardening 
Materials in Duetine Fracture Problems," AIAA Journal, Vol. 15, No. 7, 1977. 
pp. 923-931. 

14. Atluri, S. N. and Nakagaki, M. and Chen, W. H., "Fracture Analysis Under 
Large Scale Yielding : A Finite Deformation Embedded Singularity, Elastic- 
plastic Incremental Finite Element Solution". ASTM STP 631, American 
Society for Testing and Materials, 1977, pp. 43-61. 

15. Atluri, S. N., "On Some New General and Complementary Energy Theorems for 
the Rate Problems in Finite Strain, Classical Elasto-Plasticity" , Georgia 
Institute of Technology, Report GIT-ESM-SA-78-10, Aug. 1978. (In review 
for publication in J. Mechanics and Physics of Solids ). 

16. Atluri, S. N., "Complementary Energy Principles; Finite Strain Problems; 
and Finite Elements", presented at IUTAM symp. on Variational Methods in 
Engineering, Northwestern University, Sept. 11-13, 1978 (to appear in 
Proceedings to be published by Pergamon Press; S. Nemet-Nassar and K. 

Washizu, Editors). 

17. Hofmeister, L. D. , Greenbaum, G. A., and Evensen, D. A., "Large Strain, 
Elasto-Plastic Finite Element Analysis," AIAA Journal, Vol. 9, No. 7, 
pp. 1248-1254, 1971. 

18. Chitaley, A. D. and Mclintock, F. A., "Elastic-plastic Mechanics of Steady 
Crack Growth under Anti-plane Shear". Journal of Mechanics and Physics 

of Solids", Vol. 19, 1971, pp. 147-163. 

19. Rice, J. R. , "Elastic-plastic Models for Stable Crack Growth," in Mechanics 
and Mechanisms of Crack Growth (Proc . Conf . at Cambridge, England f Apr .1973) , 
ed. M. J. May, British Steel Corporation Physical Metallurgy Center Publi- 
cation, 1975, pp. 14-39. 

20. Broberg, K. B., "On Stable Crack Growth," Journal of Mechanics and Physics 
of Solids , 1975, Vol. 23, pp. 215-237. 

21. Prager, W. , "A New Method of Analyzing Stress and Strains in Work-Hardening 
Plastic Solids," Journal of Applied Mechanics , Vol. 23, 1956, p. 493. 

22. Ziegler, H. , "A Modification of Prager's Hardening Rule,", Quart. Appl . 

Math . , Vol. 17, No. 1, 1959 p. 55. 

23. Bernard, P. J., Lindley, T. C., and Richards, C. E., "Mechanisms of Over- 
load Retardation During Fatigue Crack Propagation", Fatigue Crack Growth 
Under Spectrum Loads , ASTM STP 595, American Society for Testing and 
Materials, 1976, pp. 78-97. 

24. Bowie, 0. L. and Neal, P. M., "A Note on the Central Crack in a Uniformally 
Stressed Strip", International Journal of Fracture Mechanics , Vol. 2, 

Nov. 1970, pp. 181-182. 


52 




53 



Fig. 2. 8-Noded Isoparametric "Regular" Element. 
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Fig. 3. Yield Surface Translation : Kinematic Hardening Law. 
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Fig. 5. Finite Element Model of a Center Cracked Specimen Under 
Uniaxial Cyclic Loading. ( Singular Sector Elements 
Shown Within Detail "A") 
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PORTION OF DETAIL "A" IN FIG. 5 
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Fig. 8a - 8d. Representative Patterns of Crack Closure ; (8a - c) : 
Crack Closes Only at Nodes Closest to Crack-tip ; 

8d : Crack Closure Occurs Also at Nodes Away From 
Crack-tip. 
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Fig. 9. Crack-Closure and Crack-Opening Stresses in Constant Amplitude (R=0) 
Cyclic Loading. 
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Fig. 18a - c : Crack-Closure and Crack-Opening Stress for Three 
Different Cases of a Single Over-Load. 
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Fig. 21. Crack-Line Profile During Unloading After a Single Over-Load. 
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Fig, 25. Plastic-Zone Near Crack-Tip in a Center- Cracked Panel 
Under Pure Shear Cyclic Loading, 
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Fig. 27. Maximum Opening Displacement at a Node Closest to Crack-Tip 
and Total Compressive Force over Closed Area Near Crack-Tip 
in Each Cycle, in a Single Over-Load Case. 
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